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WAVELET TRANSFORM FOR TIME-FREQUENCY 
ANALYSIS OF THE VIBRATIONAL SIGNATURE 
AND ITS APPLICATION 

By 

Jae-Jin Jeon and Young S. Shin 



ABSTRACT 

Wavelet transform is applied to the analysis of vibration 
signatures in order to verify the ability of the detection of 
abnormal condition. It can well describe the dynamics of the 
signal’s spectral composition of a non-stationary and 
stationary signal to be measured and presented in the form of 
3-D time-frequency map. Although wavelet has been developed 
over about ten years in the mathematics and physics, its 
engineering applications is a first stage. The objective of this 
report outlines the definition of the wavelet transform and is to 
discuss the properties of the wavelet transform as new tool for 
the vibration analysis, and then demonstrates how it may be 
applied to the machinery condition monitoring. 



I. INTRODUCTION 



Wavelets are a very popular topic of the signal processing and applied 
mathematics. In the last ten years, an interest in them has grown at an 
fast rate in signal and numerical analysis [Beylkin, Coifman and Rokhlin, 
1991, Heil, 1990 and Resnikoff and Burrus, 1990]. Wavelet analysis appears 
to new subject for the time-frequency analysis of the vibrational signatures. 
Traditional spectral analysis provides spectral values which are 
independent of time. It is assumed to be ergodic and stationary signal with 
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time. However the signal associated with most vibrational phenomena are 
in general time varying, which means that their characteristics change 
with time and they have various features in different time frames. For 
example, the vibration during the start-up of an engine or pump is non- 
stationary, the sound pressure generated from speaker is nonstationary, 
and so on. In case of the signal containing some transient or nonstationary 
conditions, the traditional approach in signal analysis fails to describe the 
dynamics of the signal's frequency components changes. Changes in the 
condition of a component such as a gear can be expected to cause some 
change in the vibration generated by mechanical system. Very little 
damage detection can be performed using the vibration signal directly from 
the system because the small changes generated by early damage may be 
masked by the normal vibration of the system. 

Two methods of signal analysis for nonstationary application are 
commonly used in time-frequency domain. The windowed Fourier 
transform, otherwise so called short time Fourier transform(STFT), has a 
short time window of a fixed size centered at time t as figure 1(a). The 
windowed Fourier transform is given as follows 

F(^o,t) = |g(T - t)e-'“'s(T)dT (1) 

where F(o),t) is the windowed Fourier transform, g(t) window function and 
s(t) time signal. The range of integration is from -«> to <». If the length of 
the window is time duration T, the resolution of time frequency domain 
depends on T. Its frequency bandwidth or frequency resolution is 
approximately 1/T. Therefore this method have the limitation of resolution 
in both time and frequency domain simultaneously. 

The second method, often called Wigner distribution method, is based 
on the instantaneous power spectra defined as following equation (2) 

= — ?5(/- t/ 2 ) 5(r-hT/2) dr ( 2 ) 

2n ^ 
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where >v((y,t) is Wigner distribution function. Equation (2) is called the 
Wigner distribution of 5(f)in figure 1(b). We discussed well at references 
[Jeon and Shin (1993) and Shin, Jeon and Spooner (1993)] about the 
characteristics of Wigner-Ville distribution. For a nostationary signal 
analysis, spectrogram is commonly used, which is based on the 
assumption that it is a collection of a short duration stationary signal. 











Time (sec) 

(b) Wigner distribution 

Figure 1. Calculation of the windowed Fourier transform and 
instantaneous spectral density(Wigner distribution) of a vibration signal. 
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A major drawback of this approach is that the frequency resolution is 
directly affected by the duration of short stationary time, which 
subsequently determines the time resolution. The frequency and time 
resolution of the Wigner distribution are not determined by the short 
duration but rather determined by the selection of desired resolution of the 
signal itself, but may not be appropriate for signals containing patterns at 
both large and very small scales. 

The windowed Fourier transform is well portray the characteristics of 
signals in which all of the patterns appear at approximately the same 
scale, but may not be appropriate for signals containing patterns at both 
large and small scales because of the time fixed window size. The 
multifrequency channel decomposition, which are an intermediate between 
a time and a Fourier representation, have foimd many applications in 
signal and image processing [Mallat, 1989]. Much recent research has 
been focused on this domain with the modeling of a new decomposition 
called the wavelet transform, which is presented signal by summation of 
family of functions which are the dilations and translations of a unique 
function called a wavelet. Instead of portraying a signal into harmonic 
functions ie'‘^ in Fourier transform), the signal is presented into a series of 
orthogonal basis functions of finite length. Each wavelet is located at a 
different position on the time axis. At the finest scale, wavelets may be very 
short indeed; at a coarse scale, they may be very long. Alternatively very 
small disturbances in a record of machinery vibration can be easily 
characterized from a wavelet 3-D or 2-D map in which the mean-square 
value of the time record is shown over wavelet scale and position. 

One important property of a wavelet transform is its ability to 
characterize easily the local regularity of a fimction. Simply by a change of 
the scale parameter( dilation) in the wavelet transform, many scales of local 
structure can be described by a distribution in the time-scale plane. Wang 
and McFadden (1993) introduced the advantage for examining the vibration 
signal generated by a gear and D. E. Newland (1993) well investigated the 
properties of the wavelet as a new tool for the analysis of vibration records. 
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This report outlines the definition of the wavelet transform, compare 
with the advantages and the disadvantages of wavelet and pseudo Wigner- 
Ville distribution in time-frequency domain analysis, and then 
demonstrates how it may be applied to analysis of the vibration signals for 
the machinery condition monitoring. 
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II. DEFINmON OF WAVELET TRANSFORM 



To overcome the limitation of the fixed resolution of the windowed 
Fourier transform in the time and frequency domains by decomposing the 
signal s(t) into a family of functions which are the translation and dilation 
of an unique function defined the continuous wavelet transform as 



W{a,b) = ^ 



(3) 



where W{a,b) is wavelet transform, \y is an analyzing wavelet, a 
represents a time dilation, b a time translation, and bar for complex 
conjugate. The normalization factor 1/Va is perhaps most effectively 
visualized as endowing |W(a,f))| with imit of power/Hertz [Shensa, 1992]. 



We consider the space L2(R) of measurable fxmction defined on the 



real line R (R := (-oo, »)), certain weak 'admissibility' conditions are usually 
required on ^(t) [Shensa, 1992]; 




dco < 



OO 



(4) 



where ^(fij)is the Fourier transform of They ensure that the 

transformation is a bounded invertible operator in appropriate space 
[Daubechies, 1988]. If is differentiable, then it suffices that y/ be 

zero mean, that is. 



j V^(t) dt = 0 (5) 



for equation (4) to be satisfied. In particvilar, since the local average values 
of every function in L^fR) must 'decay' to zero at ±oo, the sinusoidal(wave) 
functions (basis of Fourier transform) do not belong to L^fR). Fourier 

series representation of any functions is in L2(0,2 ti). In fact, if it looks for 
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wavelets that generate L2(R), these wavelets should decay to zero at ±<» ; 
and for all practical purposes, the decay should be very fast [Chui, 1992]. 

In signal processing, the significance of equation (3) is well understood 
by comparing it to the windowed Foiuier transform (or short-time Fourier 
transform); 



¥{(0,b) = Jg(t - b)e'^s(t)dt. (6) 

oo 

Thus, to obtain F(CO,b), one multiplies the signal by an appropriate window 
g(t) (such as Gaussian) centered at time b and then takes the Fourier 
transform. In mathematical forms, equation (6) is an expansion of the 
signal in terms of family of functions g(t-b)e‘^, which are generated from a 
single function g(t) through translations b in time and translations 0) in 
frequency. In contrast, the wavelet transform of equation (3) is an 
expansion in function - b) / a) generated by translation b and dilation 

a in time. Thus, the continuous wavelet transform is similar to windowed 
Fourier transform with a different window size for each frequency. The 
important facts of this is that, while the basis fimctions in equation (6) have 
the same time and frequency resolution at all points of the transform plane, 
those of wavelet transform have the time resolution which decrease with 
increasing a and the frequency resolution which increase with decreasing 
a width adapted to their frequency components: at high frequency \j/ are 
very narrow, while at low frequency ij/ are much broader. As a results, the 
wavelet transform is better than the windowed Fourier transform to 
analyze on very small disturbance, i.e., high frequency phenomena. This 
property can be a best advantage in signal processing since high frequency 
characteristics are generally highly localized in time whereas slowly 
varying signals require good low frequency resolution [Shensa, 1992]. 
Figure 2 shows the typical shapes of windowed Fourier transform fimctions 
and wavelet. 

The wavelet transform means that signal s(t) is characterized by 
decomposition into a set of wavelet family with series of different frequency 
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bandwidth. For the decomposition of the signal s(t) defined on real line, it 
is necessary to shift \j/ along R. Let Z denote the set of integers; 



Z = I -1, 0, 1, }. 

The simplest way for y/ to cover all of R is to consider all the integral shift 
of V'. 



V^(t - k), k € Z. 



(7) 





Figure 2. The typical shapes of (a) windowed Fourier transform function 

2 

g(t), G)) wavelet V^(0 = ^^[Daubechies, 1992]. 



8 



Next, as in Fourier transform, it must be also considered wavelets with 
different frequencies. It is considered wavelets with frequencies partitioned 
into consecutive 'octaves' (or frequency bands). For computational 
efficiency, it will be used integral power of 2 for frequency partitioning, that 
is, the small wavelets is considered as follows 

- k), j,k € Z. (8) 

From the equation (3), ^(2^t - k) is obtained from the wavelet function 
by dilation of 1/2-' and translation of k/2-’. An important particular 
case of the discrete wavelet transform which was found is that some 
wavelets exist such that ^[^ V^(2^(t orthonormal 

basis of l 2(R) [Mallat, 1989]. 

As originally proposed by Morlet et al. (1982), y/ was a modulated 
Gaussian 



V^(t) = (9) 

and this ftmction is selected to analysis the vibration of gear box for signal 
processing applications[Wang and McFadden, 1993]. The example of a 
shape for the modulated Gaussian wavelet is shown in figure 3. The 
Fourier transform of the first wavelet family y/(t/a), i.e., no translation, is 

\j/(aco) = (10) 

which has analysis frequency (OqIq. (Oq I a is a simple frequency 
parameter which determines the analyzing wavelet. We can easily see 
that equation (9) satisfies the admissibility condition equation (4). The 
analyzing bandwidth of the wavelet is proportional to l/o, thus having a 
constant relative bandwidth(BW), that is, BW/(iOQ/a) = constant. This 

feature is also reflected in the narrow time window at higher frequency(i.e., 
at smaller a ). In general, the function y/(t) is selected by its time and 

frequency localization properties [Daubechies, 1990]. 
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The scale change may be performed by substituting a =2'^ for the 
computation efficiency. In this case the wavelet family is V? V^(2^ t), 
and the definition of the wavelet transform becomes 

W(2', 6) = ^ J s(t) dt (11) 

In this report, the unique function ^(t) can be selected to be a well-behaved 
modulated Gaussian function, given by equation (9). 




Figure 3. Example of modulated Gaussiam wavelet. 
(V'CO = e'“»' 0)0 = 2k[^ ,/„ = 2 Hz) 
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III. TIME-FREQUENCY ANALYSIS 



Suppose that y/ is any basic wavelet such that both y/ and its Fourier 
transform yf are window functions which have centers and half widths of 
time and frequency domain given by / , (U , A/, Aco , respectively. Then 
the wavelet transform of an analog signal s(t) is given as follows. 

oo T' 

W{a,6) = I s(i) it (12) 

a 

— oo 

Equation (12) localizes the signal with a time window 
[b at* - aAt, b + at* + aAt], 

where the center of the window is at 6 + at and the width is given by 2aAt . 
This is called time localization in signal ginalysis. 

On the other hand, let set the Fourier transform of y/ 

rj{(o) = yf{o) + (o*), (13) 

where T\ is the Fourier transform of y/, then Tf is also a window function 
with center at (0 and width by Ad), and by Parseval identity, the wavelet 
tr£msform in equation (12) becomes [Chui, 1992] 

W(fl,^) = ^ riia(Q)-—)do) (14) 

2tc ■' a 

— oo 

where the phase shift of is determined by translation along time axis. 
Equation (14) is also localized information of spectrum s(dJ) of the signal 
s(t) with frequency window 



1 A ^ 1 A 1 

[ — - -A( 0 , — + -Aw], 

a a a a 

where the center of window is at W la and width is given by 2 Aw/ a. 

This is called frequency localization. 

From equation (12) and (14), time-frequency window of wavelet 
transform is as follows. 



ri. * AF • 1 a 1ai xa-x 

[b at - aAt, b + at + aA/]x[ — - —Aw, — + —Aw] (15) 

a a a a 

The covered area by time-frequency window is the multiplication of time 
and frequency bandwidth around the center of window. The time -frequency 
window is shown figure 4. Eventually it is considered positive frequencies, 
i.e., a > 0, the basic wavelet I// may be chosen that the center W of is a 
positive number. The ratio of the center frequency to the width of frequency 
band is given by 



(O la W* 

— = (16) 

2Aw / a 2 Aw 

which is independent of the location of the center frequency. This is called 
constant percentage band width analysis or constant-Q frequency analysis. 
In this report, the octave band is used for the analysis of vibration signal. 
The frequency window along the frequency axis narrows for large center 
frequency 0) I a and widens for small one and the time window is opposite 
to frequency window. The area of the window is a constant, given by 
4A/AW. 



1 2 



n 



0 ) 

<72 



(0 



/>, + <7,7* l>2 + <727* / 



Figure 4. Time-frequency win<lows, <7, > 02 - 
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IV. DISCRETE WAVELET TRANSFORM 



In the continuous wavelet, the family is considered 

(17) 

a 

where /? € R, a € R+ with a and \f/ is admissible. R denotes the real 
line and R+ denotes the positive real line. It is considered with discrete 
values for a and b. For the discretization of the dilation parameter, we 
select a = Qq, where m e Z, and the dilation step Oq^ 1 is fixed. For 
convenience, it assume Gq > 1. If Oq = 1, the wavelet transform may be 

similar to windowed Fourier transform. For m = 0, it seems natural as 
well to dicretize b by taking only the integer multiples of one fixed 
where is appropriately chosen so that the y/(t - n^) cover the whole 
line. We arbitrarily fix > 0 in this report. For the different values of m, 
the width of 0 ^ times the width of yf{t) defined at section 

III, so that the choice b = nb^OQ will ensure that the discretized wavelets 
at level m cover the line in the same way that ^(t - nZj^)) do. Thus it is 
selected ,a = Oq, b = nbQo!^, where m,n range over Z, and Oq > 1, bQ> 0 
are fixed; the appropriate choices for Oq and b^ depend on the wavelet y/ 
and the characteristics of signal. The discrete form of equation (17) is 
given as following equation [Daubechies, 1992] 



V^m.n 




A “ nboQQ . 
^0 



= - nl^) 



(18) 



For the computation efficiency, we assume that a = I"', that is, Gq = 2, 
where m is termed the octave of the transform. This means that the 



frequency resolution of wavelet has a octave band. The integral equation (3) 
yields a wavelet series as following equation by using equation (18). 



W(2"’, n^) - J 2m--)s(t) dt 



(19) 



At the discrete wavelet transform, the finite energy for the wavelet 
transform is not equivalent to finite energy for the wavelet series. It 
depends on the sampling grid as well as the function l//'(t). In addition, it 

often take 6 to be a multiple of a . 



W(2"’, n2"’) 




- n)s(t)dt 



( 20 ) 



A logical step in applying the theory to discrete signal is to discretize the 
integral in equation (20) as follows. 

W(2"’, n2"’) = - n)s(k) (21) 

Octave m is only output every 2^ samples. In this form the resulting 
algorithm will not be translation invariant [Mallat, 1982]. The discrete 
wavelet transform is highly not invariant under translations. In practice 
one does not use an infinite number of scales, but cuts off very low and very 
high frequencies. 



V. EXAMPLES AND DISCUSSIONS 



A signature generated by machinery involves many informations about 
its operating condition. It can be obtained the information about the 
operating condition of machinery by applying the analysis tools for the 
vibration records. Wavelet transform is a new tool that is particularly 
suited for time-frequency analysis of nonstationary or stationary signals. 
There are many advantages of using wavelet transform for both steady and 
transient signals. We will discuss the performance of the wavelet 
transform by using simple example and compare with pseudo Wigner-Ville 
distribution(PWVD) in time-frequency domain analysis. 

Before showing some examples, it is necessary to discuss how best to 
describe the results. We have foimd that 3-D map and 2-D map of wavelet 
transform are a useful presentation for many applications. The square of 
the amplitude in equation (21) has the unit power/Hz. The distribution of 
the amplitude square over the individual wavelet and position can now be 
seen as figure 5. If the length of sampled data is shorter than the wavelet 
size, the distribution at lower wavelet level, i.e., low frequency, has the 
value at only one position. 

In order to describe the results of wavelet transform, we use the 3-D 2ind 
2-D(for the complicated signal) graphics. And for graphic, the results of 
that is distributed and reduced to 256(3-D) or 512(2-D) data point along the 
time axis. Frequency axis is a log scale (octave scale or wavelet level) and 
time axis is a linear scale in figures of wavelet transform. 



A . Harmonic wave with stepwise frequency changes 

Figure 6 shows (a) the pure sine wave with stepwise frequency changes 
500 Hz, 250 Hz and 100 Hz, (b) its PWVD and (c) the wavelet transform. The 
wavelet transform and PWVD well represent the time delay and the 
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frequency components of signal. The wavelet transform is a result with 
frequency partitioned into consecutive octaves for computational efficiency. 
Then the magnitude of 100 Hz is dispersed. But the wavelet transform 
clearly describes the time delay. From these figures, we can see that the 
PWVD has the higher resolution than the wavelet transform in frequency 
line. However the computation time of PWVD is extremely higher than the 
wavelet transform. Figure 7 shows the wavelet transform of the sine wave 
with 500 Hz in time from 0.085 sec to 0.17 sec. the wavelet transform well 
represents the time delay and the frequency band of the signal. 




Time 



Figure 5. The lattice of time-frequency localization of wavelet transform 

(N = Number of data points) 



Magnitude 



2.0 



1.0 
0.0 
- 1.0 
- 2.0 

0.000 0.064 0.128 0.192 0.256 

Time (sec) 

(a) 





(b) 

Fig.6 (continued) 
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(C) 

Figure 6. Time-frequency localization of PWVD and wavelet transform: (a) 
signal s(t), (b) its PWVD and (c) its wavelet transform (fs=2000 Hz, N=512). 
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Figure 7. Wavelet transform for step sine wave (f = 500 Hz); fs = 2000 Hz, 
N=512. 
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R Swept harmoDic wave 



The effect of sweep rate on PWVD was well investigated at Jeon and 
Shin (1993), eind Shin, Jeon and Spooner (1993). Figure 8. shows (a) the 
swept cosine wave with the sweep rate 32 Hz/sec, (b) its PWVD and (c) its 
wavelet transform. From this figures, we can see that PWVD is a useful 
tool for analysis of the signal with fast frequency change in time. However 
if the records of signal is longer, the computation time may be much 
needed. The result of the wavelet transform does not clearly represents the 
sweep condition but well describes the change of frequency range with time. 




(a) 



Fig.8 (continued) 
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(b) 

Fig.8 (continued) 
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Figure 8. Time-frequency localization of the PWVD and wavelet transform: 
(a) signal s(t) = cos (2ti 32 t^), (b) its PWVD and (c) its wavelet transform (fg = 
256 Hz, N = 256). 
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C. Harmonic wave with itches 



The interesting phenomena on the signal with an abnormal component 
as a fault were investigated. Figure 9 shows (a) the harmonic wave with 
glitch at a small region, (b) is PWVD and (c) its wavelet transform. It can 
be seen that both PWVD and wavelet transform of the signal figure 9(a) well 
represents the location of each glitch and its frequency components. 
FigurelO is the result of wavelet transform in case of sampling frequency 
8192 Hz. From these figures, it can be clearly seen that the wavelet 
transform is very useful tool in the analysis of the signal required higher 
time resolution. Figure 9(c) shows the glitch components more clearly than 
PWVD because the wavelet transform has the very narrow time window in 
high frequency region. This characteristic of the wavelet transform is 
useful to detect the fault or glitch although the fault is small and to monitor 
the condition on any vibrational machinery under the steady operation 
condition. Also the wavelet transform is more effective for the analysis of 
signal which the time record length is long, since the sweep along the 
frequency line is octave step . 




Fig.9 (continued) 



Htude 




(b) 

Fig.9 (continued) 
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(c) 

Figure 9. Time-frequency localization of PWVD and wavelet transform: (a) 
the signal s(t) [Jeon and Shin, 1993], (b) its PWVD and (c) its wavelet 
transform (fs = 256 Hz, N = 256), 
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Figure 10. Time-frequency localization of wavelet transform for the signal 
Fig. 9(a) (fs = 8192 Hz, N = 8192). 
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D. Harmonic wave with pulse 



Figure 11 illustrates about the signal including small changes 
generated by early damage. In practice, this signal is not given by this 
continuous expression, but by samples, and adding a 5-function is then 
approximated by adding a constant to one sample only. The example 
signal is generated by following equation. 

s(t) = sin(2x 100 t) + sin(27i 500 1) -(■ 1.5 6(t-0.049) + 1.5 5(t-0.061) (22) 

From this figure, we can see that the wavelet transform has the great 
advantage for the analysis of the signature including the small 
disturbance. PWVD does not give the exact information about time delay of 
pulse but well represents about the frequency components of main 
signatures. The wavelet transform does riot describe the magnitude of the 
main frequency components but well represents the time delay of pulse and 
is obtained the result by very short computation time. From this example, 
the wavelet transform is very useful tool to detect the fault although that is 
very small region on time axis. 




Fig. 11 (continued) 
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(b) 

Fig. 11 (continued) 
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(c) 

Figure 11. Time-frequency localization of PWVD and wavelet transform: 
(a) the signal s(t), (b) its PWVD and (c) its wavelet transform (fs = 8192 Hz, 
N = 1024) 
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E. Actual pump signal 



The signal of actual pump was measured at the steady state speed 
without and with abnormal condition. Figure 12 show the time signal 
pattern for each condition, and figures 13 and 14 show the results of wavelet 
transform by using 3-D or 2-D graphics. Also figure 15 are the results of 
pseudo Wigner-Ville distribution. These figures show a very interesting 
results. 

From these figures, we can easily find the abnormal condition. At 13th 
frequency band of figure 13(a) and (b), the results of wavelet transform 
obviously show the different condition and pattern for each case. And also 
at 8th frequency band, the small changes of patterns can be seen. Figure 
14(a) and (b) well represent the operation condition although the shape is 
more complicate and we can see the abnormal state of pump. 

In this case, if one wavelet to have 13th frequency band is used, it gives 
good results for diagnosis and vibration condition monitoring as very short 
computation time. Figure 15 are PWVD. Also PWVD well represents the 
abnormal condition but the computation time is very long than that of 
wavelet transform. In the case of wavelet transform, it is possible to 
independently compute for each wavelet level different from PWVD. The 
difference of computation time for the wavelet transform and PWVD is 
about 100 times for the number of sample data N=2048 in the calculation of a 
whole plane. At VAX3520, the computation time is 18 seconds for wavelet 
transform and about 30 minutes for PWVD. 

Especially the wavelet transform may be very useful to detect the small 
disturbance over long record length and to analyze the signals which has a 
long time duration and intermittent abnormal condition such as non- 
rotating machinery. The advantage of the Wigner-Ville distribution is that, 
unlike the wavelet transform or the windowed Fourier transform, it does 
not introduce a reference function(such as wavelet or window fiinction in 
Eq.(3) and (1)) against which the signal has to be integrated and the short 
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time signal. The disadvantage is that the signal enters in Wigner-Ville 
distribution in a quadratic rather than linear way, which causes many 
interference phenomena as shown in references Jeon and Shin (1993), and 
Shin, Jeon and Spooner (1993). Wigner-Ville distribution may be useful in 
some application for signals which have a very short time duration. 

0.02-1 1 1 1 1 
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(b) 

Figure 12. Time patterns of actual pump signal (a) without and (b) with 
abnormal condition. 
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Fig. 13 (continued) 
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(b) 

Figure 13. Time-frequency localization 2-D map of wavelet transform for 
Fig. 12(a) and (b), respectively (fg = 10 kHz, N=2048). 
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(a) 

Fig. 14 (continued) 
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Cb) 

Figure 14. Time-frequency localization 3-D map of wavelet transform for 
Fig. 12(a) and (b), respectively (fs=10 kHz, N=2048). 
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(a) 

Fig.l5 (continued) 
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(b) 

Fig. 13 Time-frequency localization of PWVD for Fig. 12(a) and (b), 
respectively (fs = 10 kHz, N = 2048). 
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VI. CONCLUSIONS 



The wavelet transform has been investigated and applied to analyzing 
sampled signal and the actual pump signal. The results of this research 
will be valuable asset for the analysis of vibration records and condition 
monitoring of machinery. The following conclusions can be drawn: 

(1) The wavelet transform is ideally suited for portraying the wide-band 
transient or nonstationary vibration records in time-frequency domain. 

(2) The wavelet transform has a great advantage to detect the small 
disturbance of the signal. 

(3) The computation time of wavelet transform is very short in comparison 
with other time-frequency localization techniques. 

(4) The modified Gaussian wavelet was well behavior and very effective to 
analyze the vibration records. 

(5) The wavelet transform characterizes the time-frequency localization of 
the signal well and may be useful tool for the machinery condition 
monitoring. 

(6) If the wavelet level, that is, frequency band, is selected at wavelet 
transform, its results will be a useful tool for the effective pattern 
recognition of machinery diagnosis and monitoring. 
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APPENDIX. PRCXJRAMUST 

A. Program list ofWaveletTransform(FORTRAN 77) 



c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



program wavelet 



WAVELET TRANSFORM 
by 

Dept, of Mechanical Engineering 
Naval Postgraduate School 



This program is a wavelet transform by using modified 
Gaussian wavelet (octave sweep or 1/3 octave sweep). 

Variables 

n = number of input data(n must be the power of 2.) 

fs = sampling frequency 

dt = sampling time(time interval) 

bw = cutoff frequency in highpass filter 

a = dilation 

b = translation 

w = coefficient of wavelet transform I w(a,b) I 
Array 

x(i) = input data set 
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c xl(i) = filtered input data set 

c bk(i) = filter weight 

c 
c 

dimension x(8192),bk(4 100), x 1(8 192) 
complex ai,sum 
character* 16 inname, outname 
character as 
c 

ai=cmplx(0.,l.) 

pi=atan(l.)*4. 

c 

print*, '================================: 

print* 

print*,’ WAVELET TRANSFORM’ 

print* 

print*,’================================: 

print* 

print*. What is an input filename ?’ 
read(*,10) inname 

print*. What is an output filename ?’ 
read(*,10) outname 
10 format(al6) 

c 

c select the sweep method 

c 

print*, 'What do you want the sweep method ?’ 
print*,’ if 1/1 octave, input 1’ 
print*,’ if 1/3 octave, input 2’ 
read(*,*) mswp 



c 

c 

c 

c 



read the input data file 
call indata(inname,x,fs,n,nn) 
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dt=l./fs 



c 

call smean(n,x) 
c 

print*, 'Do you want to apply highpass digital' 
print*, 'filter to the original data ? (y/n)' 
read(*,20) as 
20 format(al) 

if (as.eq.'Y'.or.as.eq.'y') then 

print*,'Enter the cutoff frequency of 
print*,'the digital highpass filter (in Hz)' 
read(*,*) bw 
endif 
c 

c signal modifications 

c 

c Application of highpass digital filter 

c 

if (as.eq.'Y'.or.as.eq.'y') then 
mo=n/2 

c calculate the filter weighting 

call filter(mo,bw,dt,bk) 

write(*,*) 'finished filtering' 

c pass the highpass filter 

do 160 i=l,n 
160 xl(i)=0. 

do 200 i=l,n 

do 170 k=-mo,mo 
j=k 
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if (k.lt.O) then 
j=k*(-l) 

endif 

j=j+l 

ll=i-k 

if (ll.lt.l) then 
Il=ll+n 

elseif (ll.gt.n) then 
ll=ll-n 

endif 
bb=-bk(j) 
if (k.eq.O) then 
bb=l-bk(j) 

endif 

xl(i)=xl(i)+bb*x(ll) 

170 continue 

200 continue 

do 210 i=l,n 
210 x(i)=xl(i) 

endif 

if (mswp.eq.l) then 
nk=nn 
fac=l. 

elseif (mswp.eq.2) then 
nk=3*nn 
fac=3. 
endif 

nr=nk-int(fac)+l 

open(7,file=outname,status='new’) 

c 

c calculate wavelet transform 

c 

tini=0. 
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ttime=dt*n 
f0=2.*pi 
write(7,*) tini 
write(7,*) ttime 
write(7,*) nn, index 

do 1000 i=l,nr 

fio=2**((nk-i+l)/fac) 

a=fio*dt 

nm=int(ttime/a) 
if (nm.eq.O) nm=l 
ddt=ttime/nm 
sr=sqrt(-alog(0.01)*2)*fio 
max=int(sr) 

if (max.ge.2*n) max=2*n 

do 900 j=l,nm 

sum=cmplx(0.,0.) 
b=ddt/2.+ddt*float(j- 1 ) 
jk=int(b/dt)+l 

do 800 k=*max,max 

t=k*dt/a 
kl=jk+k 
if(kl.lt.l) then 

if (kl.lt.*n) then 
kl=kl+2*n 

else 

kl=kl+n 

endif 

endif 

if (kl.gt.n) then 

if (kl.ge.2*n) then 
kl=kl-2*n 
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else 



kl=kl-n 

endif 

endif 

sum=sum+x(kl)*cexp(-ai*fD*t)*exp(-t*t/2.)*dt 

800 continue 

w=cabs(sum)/sqrt(a) 
w=abs(w)*abs(w) 
write(7,*) ij,w 

900 continue 

1000 continue 

close(7) 
stop 
end 
c 
c 
c 

c SUBROUTINES 

c 

C 

subroutine indata(inname,x,fs,n,nn) 

dimension x(*) 
character* 16 inname 
c 

open(5,file=inname,status='old') 
read(5,*) fs,n 
do 100 j=l,n 

read(5,*) t,x(j) 

100 continue 
close(5) 
a=fs/2. 
do i=l,20 
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di=float(ii) 
if (di.ge.fl) then 
nn=i 
goto 200 

endif 
end do 

200 continue 

return 
end 
c 
c 

subroutine smean(n,x) 
c 

c This subroutine calculates and removes the mean value of 

c the sampled data, 

c 

dimension x(*) 
real meanv 
c 

asum=0. 
do 100 i=l,n 

asum=asum+x(i) 

100 continue 

meanv=asum/n 
do 200 i=l,n 

x(i)=x(i)-meanv 
200 continue 

return 

end 

c 

c 

subroutine filter(mo,b,t,bk) 
c 

c Routine generates FIR filter weights. 
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Method devised by Potter, Bickford and Glaze. 

There are a total of 2M+1 weights... filter generates M+1. 

— variables — 

t = the sampling interval in second, 
b = cutofithalf-power point) of the filter in Hz; 
must be on the range fi*om o to l/2t. 

Results are stored in bk 

- Note; in the case of highpass filter, the value of 
weight bo must use 1-bO instead of bO. 

dimension bk(*),d(3) 

data d0/0.35577019/,d(iy0.2436983/,d(2)/0.07211497/, 

♦ d(3V0.00630165/ 
pi=atan(l.)*4. 
m=mo 

first generate plain boxcar weights 

fact=2.*b*t 

bk(l)=fact 

fact=fact*pi 

do 5 i=l,m 

fi=i 

bk(i+l)=sin(fact*fi)/(pi*fi) 
trapezoidal weighting at end 
bk(m+l)=bk(m+l)/2. 

Now apply the Potter p310 window 

sumg=bk(l) 

do 15 i=l,m 

sum=dO 

fact=pi*float(i)/float(m) 
do 10 k=l,3 

sum=sum+2 . *d(k)*cos(fact*float(k )) 
bk(i+l)=bk(i+l)*sum 
sumg=sumg+2.*bk(i+l) 
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20 



ml=m+l 
do 20 i=l,ml 
bk(i )=bk(i )/sumg 
return 
end 
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B. Program List of wavelet transform(MATLAB) 



% 

% Wavelet Transform 

% 

% by 

% Dept, of Mechanical Engineering 

% Naval Postgraduate School 

% 

% Ver. 1.0 Aug. 5 1993 

% 

% By using MATLAB Ver. 4.0 

% 

% ” Variables 

% infile = 

% fs 

% n = 

% 

% avg = 

% fact = 

% x(i) = 

% z(i,j) = 

% 

load c:\users\infile 
fs=infile(l,l); 
n=infile(l,2); 
dt=l/fs; 

x=infile(2:n+l,2); 
asum=0; 

% 

% remove mean value 

% 

for jl=l:n 

asum=asum+x(j 1); 



input filename 
sampling frequency 
the number of sampled dat 
( n must be n-th power of 2.) 
mean value of sampled data 
the factor of sweep rate 
the magnitude of signal 
the result of wavelet transform 



5 1 



end 



avg=asum/n; 
for il=l:n 

x(il)=x(il)-avg; 

end 

% 

% Decide the sweep rate 
% 

fact=input(’sweep method: 1/1 octave =1, 1/3 octave = 3 ’); 

% 

a=fs/2; 
fori 1=1:20 
ii=2'^il; 
if ii >= ft 
nn=il; 
break 
end 

end 

if fact == 1 
nk=nn; 

else 

nk=fact*nn; 

end 

nr=nk-fix(fact)+l; 

% 

% calculate wavelet transform 

% 

tini=0; 

ttime=n*dt; 

fD=2*pi; 

jm=fix(ttime/(2*dt)); 

wt=zeros(nrJm); 

for il= 1 : nr 

io=2^((nk-il+l)/fact); 

a=io*dt; 

nm=fix(ttime/a); 



52 



if nm == 0 
nm = 1; 
end 

ddt=ttime/nm; 
sr=sqrt(-log(0.01)*2)*io; 
mx=fix(sr); 
if mx >= (2*n) 
mx=2*n; 
end 



for j=l:nm 
sum=0+i*0; 
b=ddt/2+ddt*(j-l); 
jk=fix(b/dt)+l; 
for k=-mx:mx 

t=k*dt/a; 

kl=jk+k; 

if kl<l 
if kl < -n 

kl=kl+2*n; 

else 

kl=kl+n; 

end 

end 

if kl > n 

if kl >= 2*n 
kl=kl-2*n; 
else 

kl=kl*n; 

end 

end 

xl=0.*i*fD*t; 

sum=sum+x(kl)*exp(xl)*exp(-t*ty2)*dt; 

end 

w=abs(sum)/sqrt(a); 

w=abs(w)*abs(w); 
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wt(il J)=ww; 



% 

% Plotting 

% 

line=128; 
2 =zeros(nr,line); 
for il=l:nr 

io=2'^((nk-i+l)/fact); 

a=io*dt; 

k=fix(ttime/a); 

ifk==0 

k=l; 

end 

if k <= line 
step=line/k; 
k2=0; 

for ii=l;k 
kl=k2+l; 
k2=fix(step*ii); 
st=k2-kl+l; 
if St <= step 
cf=step/st; 

else 

cf=st/step; 

end 



for ij=kl:k2 

z(i l,ij)=wt(i l,ii)/cfyst; 
end 
end 
else 

step=k/line; 

k2=0; 

for ii=l:line 
kl=k2+l; 
k2=fix(step*ii); 
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st=k2-kl+l; 
ssum=0; 
for ij=kl:k2 

ssum=ssum+wt(il,ij); 

end 

if st <= step 

ssum=ssum*step/st; 

else 

ssum=ssum*stystep; 

end 

z(il,ii)=ssum; 

end 

end 

end 

xma=max(z); 
xp=max(xma)*1.2; 
dt 1 =ttime/(line- 1 ); 
for k=l:line 

xval ue=(k- 1 )*dt 1 ; 
xx(k)=xvalue; 
end 

for k=l:nr 
yy(k)=k; 
end 

surf\xx,yy,z) 
xlabeK'Time (sec)') 
ylabelCFrequency step ’) 
zlabel('Amplitude') 
axis([tini ttime 1 nr 0 xp]) 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



C. Plotting Program List (3 Dimension) 

- FOTRAN 77 and CA-DISSPLA Graphic Package • 

program plot 

This program uses the graphic package CA-DISSPLA to 
plot the results of wavelet transform 



tmax 


= time record length 


tint 


= initial time 


fmin 


= start frequency step 


fmax 


= stop frequency step 


nn 


= the number of the frequency step 


fac 


= index of sweep method 


fname 


= input data filename generated by main source 
program 



dimension rr(32768),wt(50,256) 
character*25 fname 



c 

writeC*,*) 'input file name ?’ 
read(*,20) fname 
20 format(a25) 

c 

c data distribution or reduction for 3-D graphics 

c 

line=256 



open(8,file=fname,status=’old) 

read(8,*) tini 
read(8,*) tmax 
read(8,*) nn,n,fac 
nk=nn*int(fac) 
nr=nk-int(fac)+l 
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ddt=tmax/n 
do 1000 i=l,nr 

fio=2.**((nk-i+l)/fac) 

a=£io*ddt 

k=int(tmax/a) 
if(k.eq.0)k=l 
if (k.le.line) then 

step=float(line)/float(k) 

k2=0 

do 950 ii=l,k 

read(8,*) il jl,w 
kl=k2+l 
k2=int(step*ii) 
st=float(k2-kl)+l. 
if (st.le.step) then 
coe=step/st 
else 

coe=st/step 

endif 

do 930 ij=kl,k2 
wt(i,ij)=w/coe/st 
930 continue 

950 continue 

else 

step=float(k)/float(line) 

k2=0 

do 800 ii=l,line 
kl=k2+l 
k2=int(step*ii) 
sum=0. 

st=float(k2-kl)+l. 
do 750 ij=kUt2 
read(8,*) il jl,w 
sum=sum+w 
750 continue 
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800 

1000 

2000 

100 

120 

150 



if (st.le.step) then 
sum=sum*step/st 
else 

sum=sum*st/step 

endif 

wt(i,ii)=sum 

continue 

endif 

continue 

close(8) 

smax=0. 
do 2000 i=l,nr 
do 2000 j=l, line 

if (smax.le.wt(ij)) smax=wt(ij) 
continue 

write(*,*) 'maximum=’,smax 
write(*,*) 'input the maximum scale ?' 
read(*,*) fact 

do 120 i=l,nr 
do 100 j=l, line 
k=line*(i-l)+j 
rr(k)=wt(ij) 
continue 
continue 

dd=mod(nr,2) 
if (dd.eq.O) then 
nr=nr+l 

kk=(nr- 1 )*line+ 1 
kl=nr*line 
do 150 ip=kk,kl 
rr(ip)=0. 
endif 
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c 



plotting 



c 
c 

call pdev('ln03’,ieer) 
call hwshd 
call swissm 

call shdchrOO., 1,0.002,1) 
call height(0.15) 
call physord. 1,1.2) 
call area2d(5.5,6.75) 

call messag( 'WAVELET TRANSFORM $’,100,1.1,8.2) 
call blsur 

call volm3d(8.,8.,9.) 

call x3name('Time (sec) $’, 100) 

call y3name('Frequency step', 100) 

call z3name('Amplitude $’, 100) 

call zaxangOO.) 

fmax=float(nr) 

fstep=(fmax-l)/4. 

call graf3d(tini,tmax/4.,tmax,l.,fstep,fmax,0.,’SCALE',fact) 

call surmat(iT,l,line,l,nr,l) 

call end3gr(0) 

call endpl(O) 

call donepl 

end 
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D. Plottiiig Program List (2 Dimension) 

• FOTRAN 77 and CA-DISSPLA Graphic Package • 

c 

program levelplot 



c 

c This program uses the graphic package CA-DISSPLA to 
c plot the results of wavelet transform with respect to each level, 

c 

c tmax 

c tint 

c fmin 

c fmax 

c nn 

c fac 

c inname 

c 
c 

dimension x(512),y(512),wt(50,256) 
character*25 inname 



= time record length 
= initial time 
= start frequency step 
= stop frequency step 
= the number of the frequency step 
= index of sweep method 

= input data filename generated by main source 
program 



c 

write(*,*) 'input file name ?’ 
read(*,20) inname 
20 formate a25) 

c 

c data distribution or reduction for 2*D graphics 

c 

line=512 



open(8,file=fname,status='old') 

read(8,*) tini 
read(8,*) tmax 
read(8,*) nn,n,fac 
nk=nn*int(fac) 
nr=nk-int(fac)+l 
ddt=tmax/n 
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1000 i=l^ 
fio=2.**((nk*i+l)/fac) 
a=fio*ddt 
k=int(tmax/a) 
if (k.eq.O) k=l 
if (k.le.line) then 

step=float(line)/float(k) 

k2=0 

do 950 ii=l^ 
read(8,*) iljl,w 
kl=k2+l 
k2=int(step*ii) 
st=float(k2*kl)+l. 
if (st.Ie.step) then 
coe=step/st 
else 

coe=st/step 

endif 

do 930 ij=kl,k2 
wt(i,ij)=w/coe/st 
continue 
continue 
else 

step=float(k)/float(line) 

k2=0 

do 800 ii=l,line 
kl=k2+l 
k2=int(step*ii) 
sum=0. 

st=float(k2-k 1 )+ 1. 
do 750 ij=kl,k2 
read(8,*) iljl,w 
sum=sum+w 
continue 

if (st.Ie.step) then 
sum=sum*step/st 
6 1 



else 

sum=sum*stystep 

endif 

wt(i,ii)=smn 



800 


continue 




endif 


1000 


continue 




close(8) 




smax=0. 
do 2000 i=l,nr 



do 2000 j=l, line 

if (smax.le.wt(ij)) smax=wt(ij) 
2000 continue 

write(*,*) 'maximum=',smax 
write(*,*) 'input the maximum scale ?' 
read(*,*) smax 
c 

c normalizing 

c 

do 2500 i=l,nr 
do 2500 j=i,line 

2500 wt(ij)=wt(ij)/smax 

c 

c plotting 

c 

call pdev(’ln03’,ieer) 
call hwshd 
call swissm 

call shdchr(90., 1,0.002,1) 
call height(0.15) 
call page(8.5,ll.) 
call physor(1.5,1.5) 
call area2d(4. 8,4.7) 
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call xnameCTime (sec) $', 100) 

call yname( 'Frequency step', 100) 

call headinCWAVELET TRANSFORM $',100,1.1,1) 

call thkfrm(O.Ol) 

call yaxangCO.) 

call graft tini ,tmax/4. ,tmax, 1 . , 1 . ,nr+ 1 ) 
call grid(l,l) 
do 3000 i=l,nr 
do 2700 j=l,line 
x(j)=dt*float(j) 
y(j)=float(i)+wt(iJ) 

2700 continue 

call curve(x,y,line,0) 

3000 continue 

call endpKO) 
call donepl 
stop 
end 
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